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The potential energy surface of an off-lattice model protein is characterized in detail by 



the surface. The results clearly reveal the frustration exhibited by this system and explain 

why it does not fold efficiently to the global potential energy minimum. In contrast, when 

Q ' the frustration is removed by constructing a 'Go-type' model, the resulting graph exhibits 

O ' 

the characteristics expected for a folding funnel. 

1 Introduction 

The potential energy surface (PES) of an interacting system determines its structural, 
dynamic, and thermodynamic properties. Formally, the links between the PES and these 
properties are fully defined by the stationary points on the PES, its gradient (which gives 
the forces on the particles), and the partition function. However, it is only relatively 
recently that explicit connections have been sought between the overall structure of the 
PES, or potential energy 'landscape', and the behaviour of the system it describes. This 
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approach promises to provide insight into a number of fields, including protein folding, 
global optimization and glass formation. 

In the present contribution we provide a global characterization of the PES for a model 
heteropolymer, and show how this picture explains the dynamical properties observed 
in previous simulations. In the original model 'frustration' prevents efficient relaxation 
to the global potential energy minimum. However, when the frustration is removed by 
constructing the corresponding 'Go-like' model, the deep traps disappear and the result- 
ing surface resembles a funnel. The term frustration was first used in the context of 
spin glasses,'^ where it is impossible to satisfy all favourable interactions simultaneously. 
Analogous effects exist in proteins:^ a three-dimensional structure that brings together 
two mutually attractive residues may involve generating unfavourable contacts elsewhere 
('energetic frustration'), and the interconversion of two similar structures may require the 
disruption of existing favourable interactions ('geometric frustration'). 

The major difficulty in providing a fundamental explanation of structure, dynamics 
and thermodynamics in terms of the underlying potential energy surface is that the num- 
ber of stationary points grows very rapidly with the size of the system.El This growth is, 
in fact, the basis of Levinthal's 'paradox', 1^ which points out the apparent impossibility 
of a protein finding its biologically active state in a random search amongst the astro- 
nomical number of available structures. Some attempts to resolve the paradox proposed 
a reduction in the search space from the full configuration space.l^El Although it seems 
unlikely that this reduction is the solution to the paradox, there is an implicit realization 
in such approaches that, in some way, the search is not random. In terms of the energy 
landscape there are two reasons for this. Firstly, conformations have different statistical 
weights in the thermodynamic ensemble, and secondly, they are not arranged at random 
in configuration space. Levinthal's analysis assumes that the energy landscape is flat. 
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like a golf course with a single hole corresponding to the native state. ^ By constructing a 
simple model that includes an energetic bias towards the native structure, it can be shown 
that the search time on the full conformational space is dramatically reduced to physically 
meaningful scalesPE^I 

One of the first studies to consider more explicitly the organization of the energy 
landscape was that of Leopold, Montal, and Onuchicl^ These authors proposed that the 
landscape of a natural protein consists of a collection of convergent kinetic pathways that 
lead to a unique native state which is thermodynamically the most stable. Such a land- 
scape structure was termed a 'folding funnel' because it focuses the manifold misfolded 
states towards the correct target. This approach highlights the fundamental fallacy of the 
random search in Levinthal's 'paradox'. 

Funnel theory has gained widespread acceptance through its development by Wolynes 
and coworkers in terms of a free energy landscape.!^*^ The funnel can be described in 
terms of the free energy gradient towards the native structure, and the roughness — a mea- 
sure of the barrier heights between local free energy minima, which can act as kinetic 
traps. Folding is encouraged when the roughness is not large compared with the energy 
gradient. Simulations have shown that the folding ability can be measured by the ratio of 
the folding temperature, Tf, where the native state becomes thermodynamically the most 
stable, to the glass transition temperature, Tg, where the kinetics slow down dramatically 
because of the free energy barriers .I^EUrg is usually defined as the temperature at which 
the folding time passes through a certain threshold. Folding is easiest for large Tf/Tg, 
since the native state is then statistically populated at temperatures where it is kinetically 
accessible. The effect of frustration is to increase the roughness of the energy landscape 
relative to its gradient towards the native structure, thereby hindering relaxation to the 
latter. 
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We have recently showii^^'^how a new visualization of the potential energy surface 
using disconnectivity graph^reveals the features which determine relaxation of clusters 
to their global potential energy minimum. This approach has already been used by others 
to examine the energy landscape of a tetrapeptidJ^^EH and to study the effects of con- 
formational constraints in hexapeptides^^ employing an all-atom model. In the present 
contribution we analyse a coarse-grained representation of a larger polypeptide with 46 
residues. Connected sequences of minima have been reported before for this systemp^ 
and we will show how the disconnectivity graph approach provides a clearer picture of 
the relation between the energy landscape and dynamics. 

2 The Model Potential 

Intermediate in detail between lattice and all-atom models of proteins are continuum bead 
models, in which each monomer is represented by a single bead on a chain. These off- 
lattice systems have received relatively little attention in terms of landscape analysis, but 
provide a useful medium for such an approach, since atomistic representations of proteins 
are computationally demanding. 

Here we examine the effects of frustration in a model heteropolymer introduced by 
Honeycutt and Thirumalai.'^SEil These authors proposed a 'metastability hypothesis' that 
a polypeptide may adopt a variety of metastable folded conformations with similar struc- 
tural characteristics but different energies. The particular state reached in the folding 
process depends on the initial conditions. We shall see that this scenario arises from 
frustration effects intrinsic to the model, which are not expected for a 'good folder' . 

The heteropolymer has N — 46 beads linked by stiff bonds. There are three types of 
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bead: hydrophobic (B), hydrophilic (L), and neutral (N), and the sequence is 



B9N3(LB)4N3B9N3(LB)5L. 



The potential energy is given b))^ 




V = lKrY, iW - r,)^ + ^Kq £ (0, - ee)2 




+ £ £ [A,( 1 + COS (p/)+ 5,(1+ COS 3cp,)] 



i 



N-2 N \/a\ 

+"^1 E c, (^) 

,=1 j=i+2 L V';/ 



12 




) 



6 



(1) 



where rij is the separation of beads i and j. The first term represents the bonds linking 
successive beads. The bond lengths were constrained at in Ref.|2Tl but here we follow 



a and 8 are the units of length and energy defined by the last term in Eq. (HI). To put 



121 K, such as might be used for the van der Waals interactions between argon atoms. 
The second term in Eq. dU) is a sum over the bond angles, 0/, defined by the triplets of 
atomic positions r,- to r,_|_2» with = 20erad^^ and 0e = 105°. The third term is a sum 
over the dihedral angles, cp,, defined by the quartets r,; to r,+3. If the quartet involves 
no more than one N monomer then At — Bi — 1.2, generating a preference for the trans 
conformation (cp, = 180°), whereas if two or three N monomers are involved then A, = 
and Bi = 0.2. This choice makes the three neutral segments of the chain flexible and likely 
to accommodate turns. The last term in Eq. (HJ) represents the non-bonded interactions, 
and a is set equal to r^. The coefficients for the various combinations of monomer types 



Berry et al.^by replacing these constraints with stiff springs: Kr = 231.2EO ^, where 



the energy parameter in a physical context, the value of £ suggested by Berry et al!^ i 
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are as follows. 

ijeB Cij = l 

/GL, jGL,B Cij = l 
/GN, 7eL,B,N Cij = l 

with Cij = Cji and Dij = Dp. Hence, hydrophobic monomers experience a mutual van 
der Waals attraction, and all other combinations are purely repulsive, with interactions 
involving a hydrophilic monomer being of longer range. 

The global minimum of this system, which we call the BLN model, is a four-stranded 
P-barrel,'^ illustrated in Figure [T] The hydrophobic segments congregate at the core, and 
there are turns at the neutral segments. By cutting the sequence at these regions, Vekhter 
and Berry have also used this model to study the self-assembly of the ^-barrel from the 
separated strands.^S 

3 Characterizing the Energy Landscape 

The most important points on a PES are the minima and the transition states that connect 
them. A transition state is a stationary point at which the Hessian matrix has exactly one 
negative eigenvalue whose eigenvector corresponds to the reaction coordinate. Minima 
linked by higher-index saddles (the index being the number of negative Hessian eigen- 
values) must also be linked by one or more true transition states of lower energy.E3 The 
pattern of stationary points and their connectivities define the topology of the PES. 



D^j = 1 
A7 = -l 
Du = 0, 
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3.1 Exploring the Landscape 

AH the transition states in the present work were located by eigenvector-followingpffl22l 
where the energy is maximized along one direction and simultaneously minimized in all 
the others. Details of our implementation have been given before. The minima con- 
nected to a given transition state are defined by the end points of the two steepest-descent 
paths commencing parallel and antiparallel to the transition vector (i.e., the Hessian eigen- 
vector whose eigenvalue is negative) at the transition state. Rather than steepest-descent 
minimization, we have employed a conjugate-gradient method (using only first deriva- 
tives of the potential) to calculate the pathways. This technique gives similar results, and 
has the advantage of being much faster. However, it is possible for conjugate-gradient 
minimization to converge to a stationary point of higher index than a minimum. To guard 
against this problem, each optimization was followed by reoptimization with eigenvector- 
following to a local minimum. In the majority of cases, the reoptimisation converged in 
a few steps, indicating that the conjugate-gradient method had indeed found a true mini- 
mum. 

A number of similar approaches have been developed for systematically exploring 
a PES by hopping between potential wells,'^^^^! and these can be adapted to obtain a 
topographical database in several ways. Here we want to explore the energy landscape 
thoroughly, working from the global minimum upwards. In our scheme, we commenced 
at the lowest-energy known minimum, and performed an eigenvector-following search 
for a transition state along the eigenvector with the smallest non-zero eigenvalue. Having 
located a transition state, the connected minima were found by evaluating the path as 
described above. The process was then repeated, always starting at the lowest-energy 
minimum found so far, and searching along eigenvectors in both directions in order of 
increasing eigenvalue. To enable the search to explore away from the starting minimum, 
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an upper limit, n^^, was imposed on the number of eigenvectors to be searched from each 
minimum. When n^y eigenvectors had been exhausted, the search moved onto the next- 
lowest energy minimum. We note that, even if Hev is set to its maximum value of 3N — 6, 
there is no guarantee of finding all the transition states connected to a given minimum. 

The low-energy regions of the BLN model energy landscape were explored using 
Rev = 10 until 250 minima had been found. Because of the harmonic bond potential, 
following normal modes uphill in energy did not always lead to a transition state in a 
reasonable number of iterations in this system. To compensate for this problem, the value 
of Rev was raised to 20 and the search continued until a total of 500 minima had been 
found. The final number of transition states was 636. 

3.2 Visualization 

A useful visual representation of a PES is provided by the disconnectivity graph of Becker 
and Karplus.l^This technique was first introduced to interpret a structural database of the 
tetrapeptide isobutyryl-(ala)3-NH-methyl, produced by Czerminski and Elber,'^and was 
subsequently applied to study the effects of conformational constraints in hexapeptides.^ 
The method is formally expressecP^ in the language of graph theory, but can easily be 
summarized as follows. 

At a given total energy, E, the minima can be grouped into disjoint sets, called 'super- 
basins', whose members are mutually accessible at that energy. In other words, each 
pair of minima in a super-basin are connected directly or through other minima by a path 
whose energy never exceeds E, but would require more energy to reach a minimum in 
another super-basin. At low energy there is just one super-basin — that containing the 
global minimum. At successively higher energies, more super-basins come into play as 
new minima are reached. At still higher energies, the super-basins coalesce as higher 
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barriers are overcome, until finally there is just one containing all the minima (provided 
there are no infinite barriers). 

A disconnectivity graph is constructed by performing the super-basin analysis at a se- 
ries of energies, plotted on a vertical scale. At each energy, a super-basin is represented 
by a node, with lines joining nodes in one level to their daughter nodes in the level be- 
low. The choice of the energy levels is important; too wide a spacing and no topological 
information is left, whilst too close a spacing produces a vertex for every transition state 
and hides the overall structure of the landscape. The horizontal position of the nodes is 
arbitrary, and can be chosen for clarity. In the resulting graph, all branches terminate at 
local minima, while all minima connected directly or indirectly to a node are mutually 
accessible at the energy of that node. 

Visualization of the PES in terms of connectivity patterns between minima represents 
a mapping from the full configuration space onto the underlying 'inherent structures '.ISl 
Although this approach discards information about the volume of phase space associated 
with each minimum, the density of minima with energy can provide a qualitative impres- 
sion of the volumes associated with the various regions of the energy landscape. 

Some example schematic potential energy surfaces and the corresponding discon- 
nectivity graphs are illustrated in Figure |2l The first two examples demonstrate that a 
funnel-shaped landscape produces a disconnectivity graph with a single stem leading to 
the global minimum, from which branches sprout corresponding to local minima that are 
progressively cut off as the energy descends below the barriers. The contrasting nature 
of the funnels in Figures |2ta) and (b) is immediately discernible from the corresponding 
graphs, where we see that the higher barriers and lower potential energy gradient towards 
the global minimum in (a) produce long dangling branches in the disconnectivity graph. 
Figure Etc) is qualitatively different. The PES possesses a hierarchical arrangement of 
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barriers, giving rise to multiple sub-branching in the graph. The strength of the discon- 
nectivity graph in representing the topology of the PES is that it is independent of the 
dimensionality of the system, whereas schematic plots of the potential energy itself are 
restricted to one or two dimensions. 

The disconnectivity graph for the low-energy regions of the BLN model landscape is 
shown in Figure |3l using the sample of 500 minima and 636 transition states obtained 
in Section 13.11 It is immediately apparent that the PES is not a single funnel. In fact, 
it is a good example of a rough energy landscape, with repeated splitting at successive 
nodes and long descending branches. A number of low-energy structures exist which are 
separated by high barriers. Even if the barriers were not so high, there would be little 
thermodynamic driving force towards the global minimum. The fact that the attractive 
forces are of relatively long range and non-specific character means that it is possible to 
construct many significantly different structures from common motifs such as the four 
strands in the global minimum. For example, some of the low-energy minima differ only 
by the relative positions of the two purely hydrophobic strands. These can register with 
each other in a number of positions, related visually by a parallel slide. However, such a 
slide would be an unlikely mechanism because all the non-bonded interactions would be 
disrupted at once. Instead, the shortest path between such structures typically proceeds 
through over ten separate rearrangements. 

Other ways in which low-energy structures are related involve a reorientation of the 
hydrophobic strands, so that the beads which are outermost and those that come into con- 
tact in the core in Figure [T] are interchanged. Again, such a process involves many steps 
and a high barrier. The neutral turn regions can also adopt a number of configurations. 
The barriers between structures related in this way tend to be somewhat smaller, since the 
torsion potential is weaker in these regions. 
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The same structural database that is used to construct the disconnectivity graph can 
also be analysed in terms of 'monotonic sequences' of connected minima in which the 
potential energy decreases with every step P^^^ The collection of sequences leading to a 
particular minimum define what we will call a monotonic sequence basin (MSB). Whilst 
the super-basin of the disconnectivity graph is defined at a specified energy, a monotonic 
sequence basin is a fixed feature of the landscape. 

Berry et al.^ have characterized some monotonic sequences leading to the global 
minimum of the BLN model. The sequences are connected by barriers that are relatively 
low compared with the energy gradient along the sequence, leading these authors to place 
the BLN model into the category of 'structure seekers'. We note, however, that only 67 
of our sample of 500 minima lie on monotonic sequences to the global minimum, so that 
such sequences are not representative of paths to the global minimum. Furthermore, other 
low-energy minima also lie at the bottom of separate monotonic sequences of comparable 
or even larger sets of minima. Hence, this system 'seeks' only a general |3-barrel structure; 
consideration of the arrangement of the monotonic sequences shows that significantly 
different low-energy minima will be reached from different starting configurations, and 
interconversion of these minima will be relatively slow with little preference for any given 
one. 

3.3 The Effects of Frustration 

The folding characteristics of the BLN model have recently been questioned in other 
studies. Guo and Brook^used MD simulations and a histogram method to study the 
thermodynamics of folding. They identified a collapse transition to compact states with a 
peak in the specific heat, and a folding transition in terms of a similarity parameter with 
the global minimum. The free energy surface as a function of this parameter and the com- 
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pactness showed that collapse occurs before any appreciable native structure is attained, 
rather than the cooperative collapse and structuring expected for a good folder. Nymeyer 
et al|22l inferred the roughness of the energy landscape from the model's thermodynamic 
and dynamic behaviour.'*- To demonstrate the effects of frustration, they compared their 
simulations of the BLN model with a modified version in which the frustration is largely 
eliminated. We now characterize the energy landscape of this modified model. 

To remove the effects of frustration in the BLN model, all attractive interactions be- 
tween pairs of monomers that are not in contact in the native state (global minimum) are 
removed. This transformation is equivalent to setting Dij = in Eq. ([T]) for non-bonded 
pairs of hydrophobic monomers which are separated by more than 1.167 a in the global 
minimum. This change increases the heterogeneity of the interactions, since it makes 
the attractive forces more specific. The modified potential was termed 'Go-like', follow- 
ing Go and collaborators, who constructed model lattice proteins by defining attractive 
interactions between neighbouring non-bonded monomers in an assumed ground state 
structure.!^ 

Performing a survey of the energy landscape of the Go-like model as for the BLN 
model above produced 805 transition states linking the 500 low-lying minima. The dis- 
connectivity graph is shown in Figure |4| The appearance is much more funnel-like, with 
no low-energy minima separated from the global minimum by large barriers. Relaxing 
the BLN global minimum with the Go-like potential actually produces the second-lowest 
energy structure; a similar structure differing in the orientation of one of the turns lies 
slightly lower. The energy range of the disconnectivity graph is a much larger proportion 
of the global minimum well depth than in the analogous graph for the BLN model (Figure 
|3l). This range reflects the lower density of minima per unit energy in the Go-like system 
that results from the specificity of the attractive forces. The highest-energy minima in the 
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BLN sample were still relatively compact, whereas those for the Go-like model showed 
considerable unfolding of the ^-barrel. 

The plots of energy versus shortest integrated path length to the global minimum in 
Figure |5] display the difference between the BLN and Go-like energy landscapes clearly. 
For the BLN model there is little correlation between distance and energy, whereas for the 
Go-like model the energy rises with distance, as one would expect in a funnel-like land- 
scape.!^ number of individual rearrangements along the shortest paths to the global 
minimum is shown for both models in Figure |6l The distribution for the BLN model is 
broader, with some minima lying as far as 24 steps from the global minimum, in contrast 
with a maximum of 15 for the Go-like model. This reveals the greater organization of the 
Go-like energy landscape into pathways converging at the global minimum. 

A funnel-like interpretation for the Go-like model is also encouraged by the changes 
in the average properties of the individual paths between minima, as demonstrated in 
Table [T] Uphill barriers are, on average, higher and downhill barriers lower for the Go- 
like model, producing a steeper downhill gradient between minima. However, the funnel 
of the Go-like model is far from ideal. A monotonic sequence analysis shows that only 
124 of the 500 minima lie in the primary MSB, so that the relaxation from an arbitrary 
structure to the global minimum is likely to involve a number of uphill steps. 

In simulations, Nymeyer et al.l^ found that the collapse from unfolded states and the 
formation of native structure occurred cooperatively for the Go-like model, producing 
a single narrow peak in the heat capacity. They also showed that glassy dynamics, as 
measured by non-exponential relaxation from unfolded states, starts at temperatures just 
below the collapse for the BLN model, hindering the search for the native structure. In 
the Go-like model, in contrast, glassy dynamics only set in below the folding tempera- 
ture, where the global minimum still has a large equilibrium probability. These results 
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are entirely in accord with those expected from the direct characterization of the energy 
landscape presented here. 

4 Conclusions 

The disconnectivity graph analysis of the 3-colour, 46-bead model polypeptide reveals a 
frustrated energy landscape with a number of low-lying P-barrel structures in competition 
with the global potential energy minimum. Although relaxation to one of these P-barrel 
minima may be quite efficient, much longer time scales are needed for the system to 
reliably locate the global minimum, in agreement with previous simulations. 

In contrast, when the frustration is removed by changing the potential to a Go -type 
model, the landscape is transformed to one where the global minimum should be located 
easily. The competitive low-lying minima disappear following the transformation, and the 
metastable minima are organised with an energy gradient towards the global minimum. 
Our results illustrate the utility of the disconnectivity graph approach as a tool to ratio- 
nalize and predict structural, dynamic and thermodynamic behaviour from the potential 
energy surface. 
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Tables 



Table 1: Properties of individual pathways for the BLN and Go-like models, b^^ is 
the larger (uphill) barrier height between the two minima connected by transition state 
i, and is the smaller (downhill) barrier. AEf°^ = b^^ — is the energy differ- 
ence between the two minima. The angle brackets indicate averaging over the sample of 
pathways. The units of energy are e. 



Model 


BLN 


Go-like 


(b^^^P 


2.59 
0.862 

1.73 


3.07 
0.635 

2.43 
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Figures 




Figure 1: Side and end views of the global minimum of the BLN model. Hydrophobic, 
hydrophilic, and neutral beads are shaded dark grey, white and light grey, respectively. 
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Figure 2: Schematic examples of potential energy surfaces (potential energy as a function 
of some generalized coordinate) and the corresponding disconnectivity graphs. In each 
case, the dotted lines indicate the energy levels at which the super-basin analysis has been 
made, (a) A gently sloping funnel with high barriers, (b) a steeper funnel with lower 
barriers, and (c) a 'rough' landscape. 
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Figure 4: Disconnectivity graph for the Go-like model, based on a sample of 500 minima 
and 805 transition states. The energy is in units of the parameter e. 
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Figure 5: Energy of minima as a function of the integrated path length along the shortest 
path to the global minimum. Upper panel: the BLN model; lower panel: the Go-like 
model. 
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Figure 6: Distribution of the number of rearrangements along the shortest path from a 
given minimum to the global minimum for the BLN model (black) and the Go-like model 
(grey). 
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